Simulation of sound change

Functions

produce_formants <- function(target, bias = 0) {
    noise <- sample(seq(-0.2, 0.2, by = 0.05), 1)
    outcome <- target + noise + bias
    return(outcome)
}

create_lexicon <- function(lexicon_size, vowels, group_size, formants) {
    lexicon <- tibble(
        word = 1:lexicon_size,
        vowel = rep(vowels, len = lexicon_size),
        class = sample(rep(1:(lexicon_size/group_size), each = group_size,
                           len = lexicon_size)),
        frequency = rep(1:lexicon_size, each = length(vowels), len = lexicon_size),
        f1 = rep(as.list(formants), len = lexicon_size)
    )
    
    return(lexicon)
}

create_lexicon_2 <- function(lexicon_size, vowels, group_size, formants_1, formants_2) {
    lexicon <- tibble(
        word = 1:lexicon_size,
        vowel = rep(vowels, len = lexicon_size),
        class = sample(rep(1:(lexicon_size/group_size), each = group_size,
                           len = lexicon_size)),
        frequency = rep(1:lexicon_size, each = length(vowels), len = lexicon_size),
        f1 = rep(as.list(formants_1), len = lexicon_size),
        f2 = rep(as.list(formants_2), len = lexicon_size)
    )
    
    return(lexicon)
}

resample <- function(x) {
    if (length(x) == 1) {
        return(x)
    } else {
        sample(x, 1)
    }
}

populate_lexicon <- function(lexicon) {
    lexicon_size <- nrow(lexicon)
    lexicon_frequency <- lexicon$frequency
    lexicon_f1 <- lexicon$f1
    
    for (i in 1:50000) {
        word_id <- sample(1:lexicon_size, 1, prob = lexicon_frequency)
        target <- resample(unlist(lexicon_f1[[word_id]]))
        outcome <- produce_formants(target)
        lexicon_f1[[word_id]][[length(lexicon_f1[[word_id]]) + 1]] <- outcome
    }
    
    lexicon$f1 <- lexicon_f1
    return(lexicon)
}

populate_lexicon_2 <- function(lexicon) {
    lexicon_size <- nrow(lexicon)
    lexicon_frequency <- lexicon$frequency
    lexicon_f1 <- lexicon$f1
    lexicon_f2 <- lexicon$f2
    
    for (i in 1:50000) {
        word_id <- sample(1:lexicon_size, 1, prob = lexicon_frequency)
        target_1 <- resample(unlist(lexicon_f1[[word_id]]))
        outcome_1 <- produce_formants(target_1)
        lexicon_f1[[word_id]][[length(lexicon_f1[[word_id]]) + 1]] <- outcome_1
        target_2 <- resample(unlist(lexicon_f2[[word_id]]))
        outcome_2 <- produce_formants(target_2)
        lexicon_f2[[word_id]][[length(lexicon_f2[[word_id]]) + 1]] <- outcome_2
    }
    
    lexicon$f1 <- lexicon_f1
    lexicon$f2 <- lexicon_f2
    return(lexicon)
}

get_encode <- function(condition) {
    if (condition) {
        encoding_prob <- lexicon_frequency[word_id] /
            max(lexicon_frequency)
        
        encode <- sample(c(TRUE, FALSE), 1,
            prob = c(
                encoding_prob, 1 - encoding_prob
            )
        )
    } else {
        encode <- TRUE
    }
    
    return(encode)
}

sound_shift <- function(lexicon_size, vowels, group_size, formants, biased_vowel, bias, iterations, save_freq) {
    lexicon <- create_lexicon(lexicon_size, vowels, group_size, formants) %>%
        populate_lexicon()
    
    lexicon_size <- nrow(lexicon)
    lexicon_vowel <- lexicon[["vowel"]]
    lexicon_class <- lexicon[["class"]]
    lexicon_frequency <- lexicon[["frequency"]]
    lexicon_f1 <- lexicon[["f1"]]
    
    new_lexicon_f1 <- tibble(init = 1:lexicon_size)
    
    environment(get_encode) <- environment()
    
    for (iteration in 1:iterations) {
        word_id <- sample(1:lexicon_size, 1, prob = lexicon_frequency)
        
        vowel <- lexicon_vowel[word_id]
        word_class <- lexicon_class[word_id]
        
        if (vowel == biased_vowel) {
            current_bias <- bias
        } else {
            current_bias <- 0
        }
        
        #### Produce the chosen word ####
        
        target <- resample(unlist(lexicon_f1[word_id]))
        outcome <- produce_formants(target, current_bias)
        
        if (vowel == vowels[1]) {
            pool_max <- suppressWarnings(max(unlist(
                lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[2])]
            )))
            
            encode <- get_encode(outcome <= pool_max)
            
        } else if (vowel == vowels[2]) {
            if (length(vowels) == 2) {
                pool_min <- suppressWarnings(min(unlist(
                lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[1])]
                    )
                ))
                
                encode <- get_encode(outcome >= pool_min)
                
            } else {# if length(vowels) == 3
                pool_min <- suppressWarnings(min(unlist(
                    lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[1])]
                    )
                ))
                
                pool_max <- suppressWarnings(max(unlist(
                    lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[3])]
                    )
                ))
                
                encode <- get_encode(outcome >= pool_min || outcome <= pool_max)
            }
        } else {# if vowel == vowels[3]
            pool_min <- suppressWarnings(min(unlist(
                lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[2])]
                )
            ))
            
            encode <- get_encode(outcome >= pool_min)
        }
        
        #### Encode ####
        
        if (encode) {
            lexicon_f1[[word_id]][[length(lexicon_f1[[word_id]]) + 1]] <- outcome
        }
    
    
        if (iteration %% save_freq == 0) {
            column <- as.character(iteration)
            
            new_lexicon_f1 <- mutate(
                new_lexicon_f1,
                !!column := lexicon_f1
            )
            
        }
    }
    
    lexicon <- cbind(lexicon, new_lexicon_f1) %>%
        rename(`0` = f1) %>%
        select(-init) %>%
        gather(time, f1, matches("\\d")) %>%
        mutate(time = as.integer(time)) %>%
        group_by(time, word) %>%
        mutate(f1_mean = mean(unlist(f1))) %>%
        ungroup()

    
    return(lexicon)
}

sound_shift_2 <- function(lexicon_size, vowels, group_size, formants_1, formants_2, biased_vowel, bias, iterations, save_freq) {
    lexicon <- create_lexicon_2(lexicon_size, vowels, group_size, formants_1, formants_2) %>%
        populate_lexicon_2()
    
    lexicon_size <- nrow(lexicon)
    lexicon_vowel <- lexicon[["vowel"]]
    lexicon_class <- lexicon[["class"]]
    lexicon_frequency <- lexicon[["frequency"]]
    lexicon_f1 <- lexicon[["f1"]]
    lexicon_f2 <- lexicon[["f2"]]
    
    new_lexicon_f1 <- tibble(init_1 = 1:lexicon_size)
    new_lexicon_f2 <- tibble(init_2 = 1:lexicon_size)
    
    environment(get_encode) <- environment()
    
    for (iteration in 1:iterations) {
        word_id <- sample(1:lexicon_size, 1, prob = lexicon_frequency)
        
        vowel <- lexicon_vowel[word_id]
        word_class <- lexicon_class[word_id]
        
        if (vowel == biased_vowel) {
            current_bias <- bias
        } else {
            current_bias <- 0
        }
        
        #### Produce the chosen word ####
        
        target <- resample(unlist(lexicon_f1[word_id]))
        outcome <- produce_formants(target, current_bias)
        target_2 <- resample(unlist(lexicon_f2[word_id]))
        outcome_2 <- produce_formants(target_2)
        
        if (vowel == vowels[1]) {
            pool_max <- suppressWarnings(max(unlist(
                lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[2])]
            )))
            
            pool_min <- suppressWarnings(min(unlist(
                lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[2])]
            )))
            
            pool_max_2 <- suppressWarnings(max(unlist(
                lexicon_f2[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[2])]
            )))
            
            pool_min_2 <- suppressWarnings(min(unlist(
                lexicon_f2[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[2])]
            )))
            
            encode <- get_encode(
                (outcome >= pool_min & outcome <= pool_max) &
                    (outcome_2 >= pool_min_2 & outcome_2 <= pool_max_2)
            )
            
        } else if (vowel == vowels[2]) {
            if (length(vowels) == 2) {
                pool_max <- suppressWarnings(max(unlist(
                    lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[1])]
                )))
                
                pool_min <- suppressWarnings(min(unlist(
                    lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[1])]
                )))
                
                pool_max_2 <- suppressWarnings(max(unlist(
                    lexicon_f2[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[1])]
                )))
                
                pool_min_2 <- suppressWarnings(min(unlist(
                    lexicon_f2[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[1])]
                )))
                
                encode <- get_encode(
                    (outcome >= pool_min & outcome <= pool_max) &
                        (outcome_2 >= pool_min_2 & outcome_2 <= pool_max_2)
                )
                
            } else {# if length(vowels) == 3
                pool_max <- suppressWarnings(max(unlist(
                    lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[1])]
                )))
                
                pool_min <- suppressWarnings(min(unlist(
                    lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[1])]
                )))
                
                pool_max_2 <- suppressWarnings(max(unlist(
                    lexicon_f2[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[1])]
                )))
                
                pool_min_2 <- suppressWarnings(min(unlist(
                    lexicon_f2[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[1])]
                )))
                
                pool_max_3 <- suppressWarnings(max(unlist(
                    lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[3])]
                )))
                
                pool_min_3 <- suppressWarnings(min(unlist(
                    lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[3])]
                )))
                
                pool_max_2_3 <- suppressWarnings(max(unlist(
                    lexicon_f2[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[3])]
                )))
                
                pool_min_2_3 <- suppressWarnings(min(unlist(
                    lexicon_f2[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[3])]
                )))
                
                encode <- get_encode(
                    (outcome >= pool_min & outcome <= pool_max) &
                        (outcome_2 >= pool_min_2 & outcome_2 <= pool_max_2) ||
                        (outcome >= pool_min_3 & outcome <= pool_max_3) &
                        (outcome_2 >= pool_min_2_3 & outcome_2 <= pool_max_2_3)
                )
            }
        } else {# if vowel == vowels[3]
            pool_max <- suppressWarnings(max(unlist(
                lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[2])]
            )))
            
            pool_min <- suppressWarnings(min(unlist(
                lexicon_f1[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[2])]
            )))
            
            pool_max_2 <- suppressWarnings(max(unlist(
                lexicon_f2[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[2])]
            )))
            
            pool_min_2 <- suppressWarnings(min(unlist(
                lexicon_f2[which(lexicon_vowel[which(lexicon_class == word_class)] == vowels[2])]
            )))
            
            encode <- get_encode(
                (outcome >= pool_min & outcome <= pool_max) &
                    (outcome_2 >= pool_min_2 & outcome_2 <= pool_max_2)
            )
        }
        
        #### Encode ####
        
        if (encode) {
            lexicon_f1[[word_id]][[length(lexicon_f1[[word_id]]) + 1]] <- outcome
            lexicon_f2[[word_id]][[length(lexicon_f2[[word_id]]) + 1]] <- outcome_2
        }
    
    
        if (iteration %% save_freq == 0) {
            column_1 <- paste("f1", as.character(iteration), sep = "_")
            column_2 <- paste("f2", as.character(iteration), sep = "_")
            
            new_lexicon_f1 <- mutate(
                new_lexicon_f1,
                !!column_1 := lexicon_f1
            )
            
            new_lexicon_f2 <- mutate(
                new_lexicon_f2,
                !!column_2 := lexicon_f2
            )
            
        }
    }
    
    lexicon <- cbind(lexicon, new_lexicon_f1, new_lexicon_f2) %>%
        rename(`f1_0` = f1, `f2_0` = f2) %>%
        select(-init_1, -init_2) %>%
        gather(time, formant, matches("f[12]_\\d")) %>%
        separate(time, c("formant_number", "time")) %>%
        spread(formant_number, formant) %>%
        mutate(time = as.integer(time)) %>%
        group_by(time, word) %>%
        mutate(f1_mean = mean(unlist(f1)), f2_mean = mean(unlist(f2))) %>%
        ungroup()

    
    return(lexicon)
}

Two vowels, F1 (1)

Simulation

set.seed(seed)
lexicon <- sound_shift(
    lexicon_size = 100,
    vowels = c("BART", "BAT"),
    group_size = 5,
    formants = c(6.5, 5.5),
    biased_vowel = "BART",
    bias = -0.3,
    iterations = 200000,
    save_freq = 10000
) %>%
    mutate(
        freq_bin = ifelse(
            frequency < max(frequency)/2,
            "low",
            "high"
        )
    )

Plotting

lexicon %>%
    ggplot(aes(time, f1_mean)) +
    geom_jitter(alpha = 0.3, size = 0.5) +
    geom_smooth(aes(colour = vowel))

## `geom_smooth()` using method = 'gam'

lexicon %>%
    ggplot(aes(time, f1_mean, colour = frequency, group = word)) +
    geom_line() +
    facet_grid(. ~ vowel)

lexicon %>%
    ggplot(aes(time, f1_mean, colour = freq_bin, group = word)) +
    geom_line(alpha = 0.5) +
    facet_grid(. ~ vowel)

lexicon %>%
    ggplot(aes(time, f1_mean, colour = freq_bin)) +
    geom_point(alpha = 0.1, size = 0.5) +
    geom_smooth(se = FALSE) +
    facet_grid(. ~ vowel)

## `geom_smooth()` using method = 'loess'

lexicon %>%
    filter(time == 0) %>%
    ggplot(aes(frequency, f1_mean)) +
    geom_point(alpha = 0.1) +
    geom_smooth(aes(colour = vowel), method = "lm") +
    ylim(5, 7)

lexicon %>%
    filter(time == 100000) %>%
    ggplot(aes(frequency, f1_mean)) +
    geom_point(alpha = 0.1) +
    geom_smooth(aes(colour = vowel), method = "lm") +
    ylim(5, 7)

lexicon %>%
    filter(time == 200000) %>%
    ggplot(aes(frequency, f1_mean)) +
    geom_point(alpha = 0.1) +
    geom_smooth(aes(colour = vowel), method = "lm") +
    ylim(5, 7)

Three vowels, F1 (2)

Simulation

set.seed(seed)
lexicon_2 <- sound_shift(
    lexicon_size = 100,
    vowels = c("BART", "BAT", "BET"),
    group_size = 5,
    formants = c(6.5, 5.5, 4.5),
    biased_vowel = "BART",
    bias = -0.3,
    iterations = 200000,
    save_freq = 10000
) %>%
    mutate(
        freq_bin = ifelse(
            frequency < max(frequency)/2,
            "low",
            "high"
        )
    )

Plotting

lexicon_2 %>%
    ggplot(aes(time, f1_mean)) +
    geom_jitter(alpha = 0.3, size = 0.5) +
    geom_smooth(aes(colour = vowel))

## `geom_smooth()` using method = 'loess'

lexicon_2 %>%
    ggplot(aes(time, f1_mean, colour = frequency, group = word)) +
    geom_line() +
    facet_grid(. ~ vowel)

lexicon_2 %>%
    ggplot(aes(time, f1_mean, colour = freq_bin, group = word)) +
    geom_line(alpha = 0.5) +
    facet_grid(. ~ vowel)

Two vowels, F1, F2 (3)

Simulation

set.seed(seed)
lexicon_3 <- sound_shift_2(
    lexicon_size = 100,
    vowels = c("BART", "BAT"),
    group_size = 5,
    formants_1 = c(6.5, 5.5),
    formants_2 = c(12.7, 13),
    biased_vowel = "BART",
    bias = -0.3,
    iterations = 200000,
    save_freq = 10000
) %>%
    mutate(
        freq_bin = ifelse(
            frequency < max(frequency)/2,
            "low",
            "high"
        )
    )

Plotting

lexicon_3 %>%
    ggplot(aes(time, f1_mean)) +
    geom_jitter(alpha = 0.3, size = 0.5) +
    geom_smooth(aes(colour = vowel))

## `geom_smooth()` using method = 'gam'

lexicon_3 %>%
    ggplot(aes(time, f1_mean, colour = frequency, group = word)) +
    geom_line() +
    facet_grid(. ~ vowel)

lexicon_3 %>%
    ggplot(aes(time, f1_mean, colour = freq_bin, group = word)) +
    geom_line(alpha = 0.5) +
    facet_grid(. ~ vowel)

lexicon_3 %>%
    ggplot(aes(time, f1_mean, colour = freq_bin)) +
    geom_point(alpha = 0.1, size = 0.5) +
    geom_smooth(se = FALSE) +
    facet_grid(. ~ vowel)

## `geom_smooth()` using method = 'loess'

lexicon_3 %>%
    filter(time == 0) %>%
    ggplot(aes(frequency, f1_mean)) +
    geom_point(alpha = 0.1) +
    geom_smooth(aes(colour = vowel), method = "lm") +
    ylim(5, 7)

lexicon_3 %>%
    filter(time == 100000) %>%
    ggplot(aes(frequency, f1_mean)) +
    geom_point(alpha = 0.1) +
    geom_smooth(aes(colour = vowel), method = "lm") +
    ylim(5, 7)

lexicon_3 %>%
    filter(time == 200000) %>%
    ggplot(aes(frequency, f1_mean)) +
    geom_point(alpha = 0.1) +
    geom_smooth(aes(colour = vowel), method = "lm") +
    ylim(5, 7)

lexicon_3 %>%
    filter(time == 0) %>%
    ggplot(aes(f1_mean, f2_mean, colour = vowel)) +
    geom_point() +
    xlim(5, 7) + ylim(12, 13.5)

lexicon_3 %>%
    filter(time == 100000) %>%
    ggplot(aes(f1_mean, f2_mean, colour = vowel)) +
    geom_point() +
    xlim(5, 7) + ylim(12, 13.5)

lexicon_3 %>%
    filter(time == 200000) %>%
    ggplot(aes(f1_mean, f2_mean, colour = vowel)) +
    geom_point() +
    xlim(5, 7) + ylim(12, 13.5)

## Warning: Removed 1 rows containing missing values (geom_point).

Three vowels, F1, F2 (4)

Simulation

set.seed(seed)
lexicon_4 <- sound_shift_2(
    lexicon_size = 100,
    vowels = c("BART", "BAT", "BET"),
    group_size = 5,
    formants_1 = c(6.5, 5.5, 4.5),
    formants_2 = c(12.7, 13, 13.3),
    biased_vowel = "BART",
    bias = -0.3,
    iterations = 200000,
    save_freq = 10000
) %>%
    mutate(
        freq_bin = ifelse(
            frequency < max(frequency)/2,
            "low",
            "high"
        )
    )

lexicon_4 %>%
    ggplot(aes(time, f1_mean, colour = freq_bin)) +
    geom_jitter(alpha = 0.2, size = 0.5) +
    geom_smooth(se = FALSE) +
    facet_grid(. ~ vowel)

## `geom_smooth()` using method = 'loess'

lexicon_4 %>%
    filter(time == 0) %>%
    ggplot(aes(f1_mean, f2_mean, colour = vowel)) +
    geom_point() +
    xlim(4, 7) + ylim(12, 14)

lexicon_4 %>%
    filter(time == 100000) %>%
    ggplot(aes(f1_mean, f2_mean, colour = vowel)) +
    geom_point() +
    xlim(4, 7) + ylim(12, 14)

lexicon_4 %>%
    filter(time == 200000) %>%
    ggplot(aes(f1_mean, f2_mean, colour = vowel)) +
    geom_point() +
    xlim(4, 7) + ylim(12, 14)

Analysis

Two vowels, F1

lexicon <- lexicon %>%
    mutate(
        vowel_ord = ordered(vowel, levels = c("BART", "BAT")),
        word_ord = as.ordered(word)
    ) %>%
    arrange(word, time) %>%
    create_event_start("word")

lexicon_gam <- bam(
    f1_mean ~
        frequency +
        vowel_ord +
        s(time, bs = "cr") +
        s(time, bs = "cr", by = frequency) +
        s(time, bs = "cr", by = vowel_ord) +
        s(frequency, bs = "cr") +
        ti(time, frequency) +
        ti(time, frequency, by = vowel_ord) +
        s(time, word_ord, bs = "fs"),
    data = lexicon,
    method = "fREML"
)

summary(lexicon_gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## f1_mean ~ frequency + vowel_ord + s(time, bs = "cr") + s(time, 
##     bs = "cr", by = frequency) + s(time, bs = "cr", by = vowel_ord) + 
##     s(frequency, bs = "cr") + ti(time, frequency) + ti(time, 
##     frequency, by = vowel_ord) + s(time, word_ord, bs = "fs")
## 
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.8451582  0.0280640  208.28   <2e-16 ***
## frequency   -0.0017286  0.0009714   -1.78   0.0754 .  
## vowel_ord.L -0.5890718  0.0195608  -30.11   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                       edf    Ref.df        F  p-value    
## s(time)                         7.363e+00 7.524e+00   13.147  < 2e-16 ***
## s(time):frequency               5.255e+00 5.823e+00   97.520  < 2e-16 ***
## s(time):vowel_ordBAT            8.658e+00 8.823e+00   66.936  < 2e-16 ***
## s(frequency)                    1.162e-04 1.165e-04    0.418    0.994    
## ti(time,frequency)              7.312e+00 7.407e+00    6.534 1.84e-07 ***
## ti(time,frequency):vowel_ordBAT 7.135e+00 7.237e+00    9.155 1.96e-11 ***
## s(time,word_ord)                6.818e+02 9.930e+02 4847.119  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 1069/1072
## R-sq.(adj) =      1   Deviance explained =  100%
## fREML = -7684.5  Scale est. = 7.2608e-06  n = 2100

fvisgam(lexicon_gam, view = c("time","frequency"),
        cond = list(vowel_ord = "BART"),
        rm.ranef = TRUE
        )

## Summary:
##  * frequency : numeric predictor; with 30 values ranging from 1.000000 to 50.000000. 
##  * vowel_ord : factor; set to the value(s): BART. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

fvisgam(lexicon_gam, view = c("time","frequency"),
        cond = list(vowel_ord = "BAT"),
        rm.ranef = TRUE
        )

## Summary:
##  * frequency : numeric predictor; with 30 values ranging from 1.000000 to 50.000000. 
##  * vowel_ord : factor; set to the value(s): BAT. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

Three vowels, F1

lexicon_2 <- lexicon_2 %>%
    mutate(
        vowel_ord = ordered(vowel, levels = c("BART", "BAT", "BET")),
        word_ord = as.ordered(word)
    ) %>%
    arrange(word, time) %>%
    create_event_start("word")

lexicon_2_gam <- bam(
    f1_mean ~
        frequency +
        vowel_ord +
        s(time, bs = "cr") +
        s(time, bs = "cr", by = frequency) +
        s(time, bs = "cr", by = vowel_ord) +
        s(frequency, bs = "cr") +
        ti(time, frequency) +
        ti(time, frequency, by = vowel_ord) +
        s(time, word_ord, bs = "fs"),
    data = lexicon_2,
    method = "fREML"
)

summary(lexicon_2_gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## f1_mean ~ frequency + vowel_ord + s(time, bs = "cr") + s(time, 
##     bs = "cr", by = frequency) + s(time, bs = "cr", by = vowel_ord) + 
##     s(frequency, bs = "cr") + ti(time, frequency) + ti(time, 
##     frequency, by = vowel_ord) + s(time, word_ord, bs = "fs")
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.434199   0.027929 194.575  < 2e-16 ***
## frequency   -0.006797   0.001431  -4.749 2.25e-06 ***
## vowel_ord.L -1.257763   0.023874 -52.683  < 2e-16 ***
## vowel_ord.Q -0.074506   0.023995  -3.105  0.00194 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                       edf    Ref.df        F p-value    
## s(time)                         7.857e+00 7.929e+00   91.547 < 2e-16 ***
## s(time):frequency               1.601e+00 1.733e+00  364.850 < 2e-16 ***
## s(time):vowel_ordBAT            8.714e+00 8.866e+00   88.733 < 2e-16 ***
## s(time):vowel_ordBET            8.675e+00 8.846e+00   78.286 < 2e-16 ***
## s(frequency)                    3.792e-05 3.803e-05    0.074 0.99867    
## ti(time,frequency)              1.000e+00 1.000e+00  105.938 < 2e-16 ***
## ti(time,frequency):vowel_ordBAT 1.000e+00 1.000e+00    6.761 0.00941 ** 
## ti(time,frequency):vowel_ordBET 1.651e+00 1.663e+00    6.107 0.01626 *  
## s(time,word_ord)                6.623e+02 9.900e+02 4927.740 < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 1095/1098
## R-sq.(adj) =      1   Deviance explained =  100%
## fREML =  -7846  Scale est. = 6.4822e-06  n = 2100

fvisgam(lexicon_2_gam, view = c("time","frequency"),
        cond = list(vowel_ord = "BART"),
        rm.ranef = TRUE
        )

## Summary:
##  * frequency : numeric predictor; with 30 values ranging from 1.000000 to 34.000000. 
##  * vowel_ord : factor; set to the value(s): BART. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

fvisgam(lexicon_2_gam, view = c("time","frequency"),
        cond = list(vowel_ord = "BAT"),
        rm.ranef = TRUE
        )

## Summary:
##  * frequency : numeric predictor; with 30 values ranging from 1.000000 to 34.000000. 
##  * vowel_ord : factor; set to the value(s): BAT. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

fvisgam(lexicon_2_gam, view = c("time","frequency"),
        cond = list(vowel_ord = "BET"),
        rm.ranef = TRUE
        )

## Summary:
##  * frequency : numeric predictor; with 30 values ranging from 1.000000 to 34.000000. 
##  * vowel_ord : factor; set to the value(s): BET. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

plot_smooth(
    lexicon_2_gam,
    view = "time",
    cond = list(frequency = 10, vowel_ord = "BAT"),
    rug = F,
    col = "red",
    rm.ranef = TRUE
)

## Summary:
##  * frequency : numeric predictor; set to the value(s): 10. 
##  * vowel_ord : factor; set to the value(s): BAT. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

plot_smooth(
    lexicon_2_gam,
    view = "time",
    cond = list(frequency = 20, vowel_ord = "BAT"),
    rug = F,
    col = "blue",
    add = T,
    rm.ranef = TRUE
)

## Summary:
##  * frequency : numeric predictor; set to the value(s): 20. 
##  * vowel_ord : factor; set to the value(s): BAT. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

plot_smooth(
    lexicon_2_gam,
    view = "time",
    cond = list(frequency = 30, vowel_ord = "BAT"),
    rug = F,
    col = "green",
    add = T,
    rm.ranef = TRUE
)

## Summary:
##  * frequency : numeric predictor; set to the value(s): 30. 
##  * vowel_ord : factor; set to the value(s): BAT. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

Two vowels, F1, F2

lexicon_3 <- lexicon_3 %>%
    mutate(
        vowel_ord = ordered(vowel, levels = c("BART", "BAT", "BET")),
        word_ord = as.ordered(word)
    ) %>%
    arrange(word, time) %>%
    create_event_start("word")

lexicon_3_gam <- bam(
    f1_mean ~
        frequency +
        vowel_ord +
        s(time, bs = "cr") +
        s(time, bs = "cr", by = frequency) +
        s(time, bs = "cr", by = vowel_ord) +
        s(frequency, bs = "cr") +
        ti(time, frequency) +
        ti(time, frequency, by = vowel_ord) +
        s(time, word_ord, bs = "fs"),
    data = lexicon_3,
    method = "fREML"
)

summary(lexicon_3_gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## f1_mean ~ frequency + vowel_ord + s(time, bs = "cr") + s(time, 
##     bs = "cr", by = frequency) + s(time, bs = "cr", by = vowel_ord) + 
##     s(frequency, bs = "cr") + ti(time, frequency) + ti(time, 
##     frequency, by = vowel_ord) + s(time, word_ord, bs = "fs")
## 
## Parametric coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.8603978  0.0257960 227.183  < 2e-16 ***
## frequency   -0.0031338  0.0008985  -3.488 0.000502 ***
## vowel_ord.L -0.5630529  0.0180089 -31.265  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                       edf    Ref.df        F  p-value    
## s(time)                         6.947e+00 7.273e+00    8.831 5.32e-11 ***
## s(time):frequency               7.789e+00 8.152e+00  121.021  < 2e-16 ***
## s(time):vowel_ordBAT            8.534e+00 8.738e+00   99.586  < 2e-16 ***
## s(frequency)                    3.215e-05 3.230e-05    0.241    0.998    
## ti(time,frequency)              1.000e+00 1.000e+00  146.535  < 2e-16 ***
## ti(time,frequency):vowel_ordBAT 4.978e+00 5.072e+00    5.354 6.34e-05 ***
## s(time,word_ord)                6.989e+02 9.930e+02 3723.126  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 1069/1072
## R-sq.(adj) =      1   Deviance explained =  100%
## fREML = -7657.5  Scale est. = 7.464e-06  n = 2100

Three vowels, F1, F2

lexicon_4 <- lexicon_4 %>%
    mutate(
        vowel_ord = ordered(vowel, levels = c("BART", "BAT", "BET")),
        word_ord = as.ordered(word)
    ) %>%
    arrange(word, time) %>%
    create_event_start("word")

lexicon_4_gam <- bam(
    f1_mean ~
        frequency +
        vowel_ord +
        s(time, bs = "cr") +
        s(time, bs = "cr", by = frequency) +
        s(time, bs = "cr", by = vowel_ord) +
        s(frequency, bs = "cr") +
        ti(time, frequency) +
        ti(time, frequency, by = vowel_ord) +
        s(time, word_ord, bs = "fs"),
    data = lexicon_4,
    method = "fREML"
)

summary(lexicon_4_gam)

## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## f1_mean ~ frequency + vowel_ord + s(time, bs = "cr") + s(time, 
##     bs = "cr", by = frequency) + s(time, bs = "cr", by = vowel_ord) + 
##     s(frequency, bs = "cr") + ti(time, frequency) + ti(time, 
##     frequency, by = vowel_ord) + s(time, word_ord, bs = "fs")
## 
## Parametric coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.449284   0.046255 117.811  < 2e-16 ***
## frequency   -0.008690   0.002599  -3.343  0.00085 ***
## vowel_ord.L -1.244137   0.022168 -56.123  < 2e-16 ***
## vowel_ord.Q -0.115654   0.022265  -5.194 2.36e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                                     edf  Ref.df        F p-value    
## s(time)                           7.808   7.882   56.477 < 2e-16 ***
## s(time):frequency                 4.464   5.050  221.389 < 2e-16 ***
## s(time):vowel_ordBAT              8.740   8.870  142.366 < 2e-16 ***
## s(time):vowel_ordBET              8.699   8.848  141.703 < 2e-16 ***
## s(frequency)                      3.832   3.837    4.560 0.00265 ** 
## ti(time,frequency)                2.936   2.966  112.574 < 2e-16 ***
## ti(time,frequency):vowel_ordBAT   2.534   2.566    4.350 0.01134 *  
## ti(time,frequency):vowel_ordBET   2.897   2.929    5.698 0.00257 ** 
## s(time,word_ord)                669.372 990.000 3862.640 < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Rank: 1095/1098
## R-sq.(adj) =      1   Deviance explained =  100%
## fREML = -7756.3  Scale est. = 6.9601e-06  n = 2100

fvisgam(lexicon_4_gam, view = c("time","frequency"),
        cond = list(vowel_ord = "BART"),
        rm.ranef = TRUE
        )

## Summary:
##  * frequency : numeric predictor; with 30 values ranging from 1.000000 to 34.000000. 
##  * vowel_ord : factor; set to the value(s): BART. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

fvisgam(lexicon_4_gam, view = c("time","frequency"),
        cond = list(vowel_ord = "BAT"),
        rm.ranef = TRUE
        )

## Summary:
##  * frequency : numeric predictor; with 30 values ranging from 1.000000 to 34.000000. 
##  * vowel_ord : factor; set to the value(s): BAT. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

fvisgam(lexicon_4_gam, view = c("time","frequency"),
        cond = list(vowel_ord = "BET"),
        rm.ranef = TRUE
        )

## Summary:
##  * frequency : numeric predictor; with 30 values ranging from 1.000000 to 34.000000. 
##  * vowel_ord : factor; set to the value(s): BET. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

plot_smooth(
    lexicon_4_gam,
    view = "time",
    cond = list(frequency = 10, vowel_ord = "BAT"),
    rug = F,
    col = "red",
    rm.ranef = TRUE
)

## Summary:
##  * frequency : numeric predictor; set to the value(s): 10. 
##  * vowel_ord : factor; set to the value(s): BAT. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

plot_smooth(
    lexicon_4_gam,
    view = "time",
    cond = list(frequency = 20, vowel_ord = "BAT"),
    rug = F,
    col = "blue",
    add = T,
    rm.ranef = TRUE
)

## Summary:
##  * frequency : numeric predictor; set to the value(s): 20. 
##  * vowel_ord : factor; set to the value(s): BAT. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
## 

plot_smooth(
    lexicon_4_gam,
    view = "time",
    cond = list(frequency = 30, vowel_ord = "BAT"),
    rug = F,
    col = "green",
    add = T,
    rm.ranef = TRUE
)

## Summary:
##  * frequency : numeric predictor; set to the value(s): 30. 
##  * vowel_ord : factor; set to the value(s): BAT. 
##  * time : numeric predictor; with 30 values ranging from 0.000000 to 200000.000000. 
##  * word_ord : factor; set to the value(s): 1. (Might be canceled as random effect, check below.) 
##  * NOTE : The following random effects columns are canceled: s(time,word_ord)
##